FAST ITERATION OF COCYLES OVER ROTATIONS AND 
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Abstract. In this paper, we develop numerical algorithms that use small requirements 
of storage and operations for the computation of hyperbolic cocycles over a rotation. 
We present fast algorithms for the iteration of the quasi-periodic cocycles and the com- 
putation of the invariant bundles, which is a preliminary step for the computation of 
invariant whiskered tori. 
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1. Introduction 

The goal of this paper is to describe efficient algorithms to compute quasi-periodic 
cocycles over rotations. We present fast algorithms for the iteration of cocycles over 
rotations and for the calculation of their invariant bundles. The main idea is to use a 
renormalization algorithm which allows to pass from a cocycle to a longer cocycle. 
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The calculation of invariant bundles for cocycles is a preliminary step for the calculation 
of whiskered invariant tori. Indeed, these algorithms require the computation of the 
projections over the linear subspaces of the linear cocycle. 

2. Some standard definitions on cocycles 

Given a matrix-valued function M : — )■ GL{2d,W) and a 

vector w G M^, we define the cocycle over the rotation T^^ associated to the matrix M by 
a function : Z x -> GL{2d, M) given by 

{M{e + {n-l)uj)---M{e) n>l, 
Id n = 0, (2.1) 

M-\6 + {n + l)u)---M-^{9) n < 1. 

Equivalently, a cocycle is defined by the relation 

Mio,e) = Id, 

M{l,0) = M{e), (2.2) 

M{n + m,e) = M{n,T^{e))M{m,e). 

We will say that M is the generator of M. Note that if M(T^) C G where G C 
GL{2d,R) is a group, then M{Z,T^') C G. 

The main example of a cocycle in this paper is 

M{e) = {DFoK){e), 

for K a parameterization of an invariant torus satisfying the invariance equation 

F o K = K oT^. 

Other examples appear in discrete Schrodinger equations jPui02j . In the above men- 
tioned examples, the cocycles lie in the symplectic group and in the unitary group, 
respectively. 

Similarly, given a matrix valued function M{6), a continuous in time cocycle Ai{t,6) 
is defined to be the unique solution of 

j^M{t,e) = M{e + oot)M{t,e), ^23) 
M{o,e) = id. 

From the uniqueness part of Cauchy-Lipschitz theorem, we have the following property 

M{e,t + s) = M{e + ujt, s)M{e,t), 

(2.4) 

Mi9,0)=ld. 
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Note that (12 .Sp and (12 .4^ are the exact analogues of (12.11) and (I2.2p in a continuous 
context. Moreover, if M(T^) C where is a sub-algebra of the Lie algebra of the Lie 
group G, then M{R, C G. 

The main example for us of a continuous in time cocycle will be 

M{e) = {DXoK){e), 

where is a solution of the invariance equation 

d^K = XoK 

and X is a Hamiltonian vector field. In this case, the cocycle A4{6,t) is symplectic. 

3. Hyperbolicity of cocycles 
One of the most crucial property of cocycles is hyperbolicity (or spectral dichotomies) 



as described in jMSMj ISSTil ISSTGaj ISS76bl [Sac78]. 

Definition 1. Given < A < /i we say that a cocycle Ai{n,6) (resp. Ai{t,6)) has a 
A, /i— dichotomy if for every 6 there exist a constant c > and a splitting depending 
on 9, 

which is characterized by: 

{xe,v) e£' ^ \M{n,d)v\ < cA"|w| , Vn > 

(3.1) 

{xe,v) eS"" ^\M{n,e)v\<c^''\v\ , Vn < 
or, in the continuous time case 

(xe, v)eS' ^ \M{t, e)v\ < cA>| , Vt > 

(3.2) 

{xe,v)e£''^\M{t,d)v\<c^'\v\, Vt < 0. 

The notation £^ and E"^ is meant to suggest that an important case is the splitting 
between stable and unstable bundles. This is the case when A < 1 < /i and the cocycle 
is said to be hyperbolic. Nevertheless, the theory developed in this section assumes only 
the existence of a spectral gap. 

In the application to the computation of tori, M{9) = {DF o K){6) and xe = K(6). 
The existence of the spectral gap means that at every point of the invariant torus K{6) 
one has a splitting so that the vectors grow with appropriate rates A, fi under iteration of 
the cocycle. In the case of the invariant torus, it can be seen that the cocycle is just the 
fundamental matrix of the variational equations so that the cocycle describes the growth 
of infinitesimal perturbations. 
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It is well known that the mappings 9 — )■ £^'^ are C"^ if M(., ) G C'" for r G N U {oo, a;} 
[HLOGcj . This result uses heavily that the cocycles are over a rotation. 

A system can have several dichotomies, but for the purposes of this paper, the definition 
[Hwill be enough, since we can perform the analysis presented here for each gap. 

One fundamental problem for subsequent applications is the computation of the in- 
variant splittings (and, of course, to ensure their existence). The computation of the 
invariant bundles is closely related to the computation of iterations of the cocycle. 

The first algorithm that comes to mind, is an analogue of the power method to compute 
leading eigenvalues of a matrix. Given a typical vector (xg, v) G we expect that, for 
n ^ 1, A4.{n,9)v will be a vector in . Even if there are issues related to the 9 

dependence, this may be practical if £^ is a 1-dimensional bundle. 

4. Equivalence of cocycles, reducibility 

Reducibihty is a very important issue in the theory of cocycles. We have the following 
definition. 

Definition 2. We say that a cocycle M{9) is equivalent to another cocycle M{9) if there 
exists a matrix valued function Q : — )■ GL(2(i, R) such that 



It is easy to check that M being equivalent to M is an equivalence relation. 
If M is equivalent to a constant cocycle (i.e. independent of 6*), we say that M is 
"reducible." 

The important point is that, when (14. ip holds, we have 



so that the iterations of reducible cocycles are very easy to compute. 

We will also see that one can alter the numerical stability properties of the iterations 
of cocycles by choosing appropriately Q. In that respect, it is also important to mention 
the concept of "quasi-reducibility" introduced by Eliasson [EliOlj . 



M{9) = Q{9 + uj)-^M{9)Q{9). 



(4.1) 



M{n, 9) = Q{9 + ncu)-^M{n, 9)Q{9). 



(4.2) 



In particular, if M is a constant matrix, we have 



Min, 9) = Q-\9 + nu)M''Q{9), 



fast numerical algorithms 
5. Algorithms for fast iteration of gogygles over rotations 
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In its simplest form, the algorithm for fast iteration of cocyles is: 



Algorithm 1 (Iteration of cocycles 1). Given M{6), compute 



M{e) = M{9 + u)M{9). 



(5.1) 



Set M — )■ M, 2uj ^ u and iterate the procedure. 

We refer to M as the renormalized cocycle and the procedure as a renormalization 
procedure. 

The important property is that applying k times the renormalization procedure de- 
scribed in Algorithm [1] amounts to compute the cocycle M.{2^,9). 

Then, if we discretize the matrix M{9) taking N points (or N Fourier modes) and 
denote by C{N) the number of operations required to perform a step of Algorithm [H we 
can compute 2^ iterates at a cost of kC{N) operations. 

Notice that the important point is that the cost of computing 2'' iterations is propor- 
tional to k. Of course, the proportionality constant depends on N. The form of this 
dependence on N depends on the details on how we compute the shift and the product 
of matrix valued functions. 

There are several alternatives to perform the transformation (15.11) . The main difficulty 
arises from the fact that, if we have points on a equally spaced grid, then 9 + u will not 
be in the same grid. We have at least three alternatives: 

(1) Keep the discretization by points in a grid and compute M{9-\-uj) by interpolating 
with nearby points. 

(2) Keep the discretization by points in a grid but note that the shift is diagonal 
in Fourier space. Of course, the multiplication of the matrix is diagonal in real 



(3) Keep the discretization in Fourier space but use the Cauchy formula for the prod- 



space. 



uct. 



The cost factor of each of these alternatives is, respectively. 



Ci{N) = 0(A), 
C2{N) = O(AlogA) 
CsiN) = 0(A2). 



(5.2) 
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Besides their cost, the above algorithms may have different stabihty and roundoff prop- 
erties. We are not aware of any study of these stability or round-off properties. The 
properties of interpolation depend on the dimension. 

In each of the cases, the main idea of the method is to precompute some blocks of the 
iteration, store them and use them in future iteration. One can clearly choose different 
strategies to group the blocks. Possibly, different methods can lead to different numerical 
(in) stabilities. At this moment, we lack a definitive theory of stability which would allow 
us to choose the blocks. 

Next, we will present an alternative consisting of using the QR decomposition for the 
iterates. As described, for instance in |Ose68t lERSSj ID WO 2] . the QR algorithm seems 
to be rather stable to compute iterates. One advantage is that, in the case of several 
gaps, it can compute all the eigenvalues in a stable way. 

Algorithm 2 (Computation of cocycles with QR). Given M{6) and a QR decomposition 

ofM{e), 

M{e) = Q{e)R{e), 

perform the following operations: 

• Compute S{e) = R{e + uj)Q{e) 

• Compute pointwise a QR decomposition of S, S{6) = Q{9)R{9). 

• Compute Q{e) = Q{e + u)Q{e) 

R{e) = R{e)R{e + io) 
me) = Q{e)R{e) 

• Set M ^ M 

R^R 
Q^Q 

and iterate the procedure. 

Since the QR decomposition is a fast algorithm, the total cost of the implementation 
depends on the issues previously discussed (see costs in (15.21) ). Instead of using QR 
decomposition, one can also perform a SVD decomposition (which is somewhat slower). 

In the case of one-dimensional maps, one can be more precise in the description of the 
method. Indeed, if the frequency u has a continued fraction expansion 

u = [ai,a2, . . . ,a„, . . .], 
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it is well known that the denominators g„ of the convergents of a; (i.e. Pn/<in = [^i, • ■ • , ^n]) 
satisfy 

Qn = (^nQn-l + Qn-2, 

qi = ai, 
qo = 1- 

As a consequence, we can consider the following algorithm for this particular case: 

Algorithm 3 (Iteration of cocycles ID). Given u = [ai, . . . , an, ■ ■ ■] and the cocycle over 
generated by M{6), define 

tj° = 0, uj^ = uj, M°(^)=Id, M^{e) = M{e + {ai-l)uj)---M{e). 

Then, for n > 2 

M(")(^) = M("~i)(^ + (a„ - ■ ■ ■ M^''-^\e)M^''-^\9) 

is a cocycle over 

and we have 

M{qn,9)=M^^Hl,9). 

The advantage of this method is that the effective rotation is decreasing to zero so 
that the cocycle is becoming in some ways similar to the iteration of a constant ma- 
trix. This method is somehow reminiscent of some algorithms that have appeared in the 
mathematical literature |Ryc92 [KriQQbl [KriQQaj . 



6. The "straddle the saddle" phenomenon and preconditioning 

The iteration of cocycles has several pitfalls compared with the iteration of matrices. 
The main complication from the numerical point of view is that the (un) stable bundle 
does depend on the base point. 

In this section we describe a geometric phenomenon that causes some instability in the 
iteration of cocycles. This instability -which is genuine- becomes catastrophic when we 
apply some of the fast iteration methods described in Section [51 The phenomenon we 
will discuss was already observed in |HL06a] . 

Since we have the inductive relation, 

M{n,e) = M{n-l,e + uj)M{e), 
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we see that we can think of computing J^{n,9) by applying J^{n — 1,9 + u) to the 
column vectors of M{6). 

The j^'^-column of M, which we will denote by mj{9), can be thought geometrically 
as an embedding from to M^'' and is given by M{6)ej where Cj is the j*'^ vector of the 
canonical basis of M^'^. If the stable space oi M.{n — 1,0 + u) has codimension ^ or less, 
there could be points 9* G such that mj{6*) G £^ ^ and such that for every 6 9* we 
have mj{e) ^ £^^. 

Clearly, the quantity 

M{n-l,e* +uj)mj{e*) 

decreases exponentially. Nevertheless, for all ^ in a neighborhood of 6* such that 6^6* 

M{n-l,e + uj)mj{9) 

will grow exponentially. The direction along which the growth takes place depends on 
the projection of mj{9) on £xe+^- 

For example, in the case d = 2, i = 1 and the stable and unstable directions are one 
dimensional, the unstable components will have different signs and the vectors M{n — 
1,^ + uj)mj{6) will align in opposite directions. An illustration of this phenomenon 
happens in Figure [H 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Figure 1. The straddle the saddle phenomenon. We plot one of the 
components of the cocycle Ai{2'',6) for the values k = 0,3,4. The case 
k = was scaled by a factor 200. 



The transversal intersection of the range of mj{6) with is indeed a true phenomenon, 
and it is a true instability. 
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Unfortunately, if mj{6) is very discontinuous as a function of 6, the discretization in 
Fourier series or the interpolation by splines will be extremely inaccurate so that the 
Algorithm [1] fails to be relevant. 

This phenomenon is easy to detect when it happens because the derivatives grow 
exponentially fast in some localized spots. 

One important case where the straddle the saddle is unavoidable is when the invariant 
bundles are non-orient able. This happens near resonances (see |HL07] ). In |HL07] . it 
is shown that, by doubling the angle the case of resonances can be studied comfortably 
because then, non-orientability is the only obstruction to the triviality of the bundle. 

6.1. Eliminating the "straddle the saddle" in the one-dimensional case. For- 
tunately, once the phenomenon is detected, it can be eliminated. The main idea is that 
one can find an equivalent cocycle which does not have the problem (or presents it in a 
smaller extent). 

In more geometric terms we observe that, even if the stable and unstable bundles are 
geometrically natural objects, the decomposition of a matrix into columns is coordinate 
dependent. Hence, if we choose some coordinate system which is reasonably close to the 
stable and unstable manifolds and we denote by Q the change of coordinates, then the 
cocycle 

Mie) = Q{e + uj)-^M{e)Q{e), 

is close to a constant. Remark that this is true only in the one- dimensional case. The 
picture is by far more involved when the bundles have higher rank. 

This may seem somewhat circular, but the circularity can be broken using continuation 
methods. Given a cocycle which is close to constant, fast iteration methods work and 
they allow us to compute the splitting. Then if we have computed Q for some M, we 
can use it to precondition the computation of neighboring M. 

The algorithms for the computation of bundles will be discussed next. 

6.2. Computation of rank-1 stable and unstable bundles using iteration of 
cocycles. The algorithms described in the previous section provide a fast way to iterate 
the cocycle. We will see that this iteration method, which is similar to the power method, 
gives the dominant eigenvalue Xmaxi^) and the corresponding eigenvector m(6'). 

The methods based on iteration rely strongly on the fact that the cocycle has one 
dominating eigenvalue which is simple. 



10 GEMMA HUGUET, RAFAEL DE LA LLAVE AND YANNICK SIRE 

Consider that we have performed k iterations of the cocycle (of course we perform 
scahngs at each step) and we have computed A4{n, 9), with n = 2^ . Then, one can easily 
read the dominant rank-1 bundle from the QR decomposition of the cocycle A^(n, 
just taking the column of Q associated to the largest value in the diagonal of the upper 
triangular matrix R. One obtains a vector m{6 + 2'^uj) (and therefore m{6) by performing 
a shift of angle —2^uj) of modulus 1 spanning the unstable manifold. Since, 

we have then 

\maM = {\M{e)m{e + u)Y[M{e)m{e + uj)]f'\ 

As it is standard in the power method, we perform scalings at each step dividing all 
the entries in the matrix M{9) by the maximum value among the entries of the matrix. 

Hence, for the simplest case that there is one dominant eigenvalue, the method pro- 
duces a section m (spanning the unstable subbundle) and a real function \maxi which 
represents the dynamics on the rank 1 unstable subbundle, such that 

M{e)m{e) = \^aMT^{o + uj). 

Following |HL06b] . under certain non-resonant conditions which are satisfied in the 
case of the stable and unstable subspaces, one can reduce the 1-dimensional cocycle 
associated to M and w to a constant. Hence, we look for a positive function p{9) and a 
constant yU such that 

\maMp{0)= ^^p{e + UJ). (6.1) 

If ^max{d) > 0, we take logarithms on both sides of the equation (16. ip . This leads to 

hgXmaxiO) +logp{9) = log IJ + log p{9 + uj), 

and taking log fj, to be the average of log Xmax{0) the problem reduces to solve for logp{9). 
The case Xmaxi^) < is analogous. Of course, p{9) and /i can be obtained just exponen- 
tiating. 
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